Many-body position operator in lattice fermionic systems with periodic boundary 

conditions 



Balazs Hetenyi 

Institut fiir Theoretische Physik, Technische Universitdt Graz, Petersgasse 16, A-8010 Graz, Austria 

and 

Mathematisches Institut, Fakultdt fiir Mathematik, Informatik und Statistik, Ludwig Maximilians Universitdt, 

Theresienstrasse 39, Miinchen 80333, Germany 
(September 1, 2009) 

A total position operator X in the position representation is derived for lattice fermionic sys- 
tems with periodic boundary conditions. The operator is shown to be Hermitian, the generator of 
translations in momentum space, and its time derivative is shown to correspond to the total current 
operator in a periodic system. The operator is such that its moments can be calculated up to any 
order. To demonstrate its utility finite size scaling is applied to the Brinkman-Rice transition as 
well as metallic and insulating Gutzwiller wavefunctions. 
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The position operator and its moments give important information about localization in quantum systems. As was 
shown by Kohn [1] metals and insulators are distinguished by the extent of their localization. Many real systems 
are periodic, and in many model systems periodic boundary conditions are imposed. In such cases the Hilbert space 
which that forms the domain of operators is restricted hence the position operator is ill-defined [2] . The single particle 
position operator in the crystal momentum representation was derived by Blount [2] and discussed extensively in the 
context of band-theory. In the crystal momentum representation this operator can be generalized to the many-body 
case [3]. To calculate the total position in the position representation Resta [4,5] suggests to average the quantity 
rpj^g expectation value of the total position operator is then defined as 

(a;) = ;^Im ln(*|e*^|*). (1) 
27r 

Via first order perturbation theory, Resta also shows [4] that the time derivative of the polarization operator based 
on the above definition gives the total current in the limit L ^ so. This idea has been applied to lattice fermionic 
systems at half-filling [5], and extended to systems at arbitrary fillings [6]. A related formalism due to Souza et 
al. [7] based on the cumulant generating function (of which Eq. (1) is a special case) establishes relations between 
localization and polarization. 

It is important to note that the position operator in this method is calculated indirectly, by first evaluating the 
expectation value of e*^. Eq. (1) is valid as can be shown [5] but the calculation of higher moments is not straight 
forward, the spread functional suggested by Resta and Sorella [5] (based on eq. (1)) is valid in the thermodynamic 
limit. 

Here it is shown that a total position operator for a lattice fermionic system with periodic boundary conditions can 
be defined as the generator of total momentum shifts. It is also demonstrated that the time derivative of the total 
position operator gives the current for a system with any number of sites (finite L). The total position operator derived 
below is such that expectation values of arbitrary powers are readily evaluated, hence an accurate assessment and 
finite size scaling of localization is enabled (up to any desired order) . The utility of the operator is then demonstrated 
via variational calculations on the Hubbard model [8-10] based on the Gutzwiller wavefunction [10,11]. 

The derivation of the total position operator is closely related to that of the total momentum operator in Ref. [12]. 
The class of models for which the formalism presented below are those used in strongly correlated systems consisting 
of site to site hopping terms and some interaction terms. An example of a lattice model is the Hubbard Hamiltonian, 

H=-tJ2 (4a(^j- + -^-c-) + "'T"n, (2) 

{i,j)(T io- 

consisting of L sites. In the following, the total position operator will be derived for the one-dimensional Hubbard 
model. Generalizations to higher dimensions and other lattice models will be discussed below. 

The real-space (Wannier state) and reciprocal-space (Bloch state) creation operators are related in the usual way, 

where Xj is the position of site j. In order to define a total position operator we first define a momentum permutation 
operator as 

P« = l- (4-4)^-20, (4) 
where c^. creates a particle in the Bloch state k. A momentum space shift operator can be defined as 

itn = Pn-ln----Pl2, (5) 



with the property that 



u X ^ J Cfe-iili,, k = 2,...,L , , 



For systems with spin- 5 particles we can define the compound momentum space shift operator as 

il = ilL^iiLl, (7) 
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with the property 

ilcj> = e'~f^Cj>ii, (8) 

where Cj^rr is an annihilation operator for particles at site Xj with spin a. 

We define the total position operator X through three conditions. First we require it to be the generator of total 
momentum shifts, i.e. 

il=e'^. (9) 

We also require X to be Hermitian, 

X = Xt. (10) 
and that the time derivative of X give the total current, 

eX = ie[H,X] = J, (11) 

which for the Hubbard model is defined as 

J = -iet {c\„Cj„ - c]„Cia). (12) 
In order to derive the explicit form of X we first define 

L-l 

g{a) = Y,ie-'"-^, (13) 

x=0 

which can be evaluated via the geometric sum formula to give 

9{a)=i- (14) 

1 — e * 

We can take the derivative of g{a) at some integer value m for a, 

g'{m) = ^Yxe-'^. (15) 

x=0 

Inverting the Fourier series, we can obtain an expression for the position x valid for a; = 0, i — 1, 

L 



For m ^ L, 



g'{m) = 27r/(e-*^ - 1), (17) 

and g'{L) can be evaluated from Eq. (15) using the arithmetic sum formula giving g'{L) = ■jt{L — 1). Thus, an overall 
expression for x reads as 



L-l 

x — 

m=l 



The right hand side of Eq. (18) is the sawtooth function /(x) = xmodL. Wc propose to take the sawtooth function as 
the definition of our position operator. Based on Eq. (9) we write the total position operator X for a many-particle 
system as a power series in the momentum shift operator as 
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It is to be emphasized that X is a genuine many-body operator (as is that of Resta [4]). 

Having defined our total position operator, we can now test whether it satisfies the requirements (Eqs. (9), (10), 
and (11)). Letting X operate on an arbitrary Wannier state (|x, cr) = c|^^ ...cjy |0)) for a system gives the result 

L-l - 2x^(xi + .-. + x,^0 \ 

Tn=l \ / 

= {{xi+ ...+XL)modL)\x,a) (20) 
where we have used Eqs. (8) and (18). Since 

H|x,a) = e' ^^|x,a), (21) 

Eq. (9) follows. Hermiticity of X follows from the unitarity of il and from the fact that il^ = 1. 

To demonstrate that the operator X satisfies the condition in Eq. (11), we first note that il commutes with the 
interaction part of the Hamiltonian. This can be shown using Eq. (8). Thus our task consists of evaluating the 
commutator [T, X] , T denoting the kinetic part of the Hubbard Hamiltonian. We first define an operator 

^=E-it-- (22) 

\ e L — 1 

The last term in the sum is divergent. However, below we show that this divergence disappears for the commutator 
[T,Y]. 

We first evaluate the commutator 

{T,Y]=f: [^Lf] (23) 
1 e * i- — 1 

We split the kinetic energy in two parts as 

= E SV'^-' (24) 



thus we can rewrite Eq. (23) as 



— ; e — 1 

Each commutator in Eq. (25) can be evaluated using Eq. (8). We obtain 

[y4,it™] = (e"'^ - l)it'"A 

[A^il™] = (1 - e-'^)A^ii"', (26) 
giving a new expression for the commutator 

L 

[T, = E - (27) 

m—l 

We now substitute the condition in Eq. (9) and we obtain 

L 

[T,Y] = E e'—^ A - A^e'—i^ . (28) 



m=l 
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It is easily seen that this commutator is zero, since X operating on a Wannier state gives an integer and 



L 

; 22I2L 



^ = 0. (29) 



m=l 



On the other hand, using the same reasoning we used to arrive at Eq. (27) it can be shown that 

L-\ 

[T, X] = ^ it^yl - Atit" , (30) 

m=l 

hence, from Eq. (27) we see that 

[T,X] = A^-A, (31) 

since il^ = 1. From Eq. (31) the expression for the current (Eq. (12)) follows straightforward. 

The total position operator X derived above can be generalized to many dimensions as follows. In higher dimensions 
the operator becomes a vector operator. The generalization of the above derivation has to be based on a generalized 
total momentum shift operator consisting of the product of all one-dimensional momentum shift operators in a 
particular direction. For example, for a three dimensional system with dimensions x,y, z a total momentum shift 
operator for the x direction (spinless case) would consist of the product of all one dimensional momentum shift 
operators 

2^JL,x=^<x^ (32) 

where ll^^^f^ denotes the total momentum shift operator in the a;-direction for a given set of coordinates y, z (Eq. (5)). 
Such an operator satisfies the commutation relation 

~ _ J Ck^-i^ky,k^^L,x, kx = ^2, L; ky, kz = 1, L , , 

^L,.Ck^,ky,k. - I CL,k, ,^WL.^^, kx = 1; ky, fc, = 1, L. 

Subsequent construction of a total position operator for a three dimensional systems follows the same steps as the 
one-dimensional case. The total momentum shift operator for a spin-^ system can be written as 

m=^L,iA^L,^,l, (34) 

where SUj is a vector operator, and i = x,y, z. A particular component of the total position operator can then be 
written as 



L-l 

m=i ^2 e- 

The commutator of operator Ri will give the current in the i direction. This is a consequence of the fact that the 
operator 211; commutes with the hoppings in directions other than i included in the Hubbard Hamiltonian. 

Extensions of the Hubbard model can also be handled. More complex interaction types (nearest neighbor, etc.) 
follow the same derivation as above, as the expression for the current does not change in this case. For more complex 
hoppings the expression for the current is modified to include the new hoppings, but the derivation presented above 
is still valid. 

For impurity models [13,14] the strategy of derivation of a total position operator is modified slightly. For example, 
the one-dimensional periodic Anderson model, in which each site contains a set of localized /-orbitals, can be written 

as 

with 

H' = Y,{Vifl,.Ci,. + H.c.}. (37) 
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In Eqs. (36) and (37) nf{i,l,a){fj ^ ^) denotes the density(creation operator) of /-orbital with label I at site i and 
with spin a. Each lattice site contains a set of / orbitals, but there are no inter-site hoppings between the localized 
/-orbitals on different sites. As a consequence the current operator is the same as that of the Hubbard model, inspite 
of the fact that the charge density includes the /-orbital terms [15]. One could construct a total position operator 
which does not include impurity orbitals, and has the same form as X derived above (only electrons in the conduction 
band enter the definition). As conduction takes place only on the standard lattice sites, not the ones associated with 
the /-orbitals, such an approach may in some cases be sufficient to characterize localization phenomena associated 
with metal-insulator transitions. However it is also possible to construct a total position operator valid for a system 
with the periodic Anderson Hamiltonian. 

To do this one has to consider the /-orbitals as separate lattices, and construct a total momentum shift operator 
for each set of /-orbitals localized on different lattice sites. One can construct an operator 



^f=QtiL--Q^il (38) 



where 

■,(0 



Qjk='^-ifli-fk,i)ihi-fk,i)- (39) 
fj^i denotes the Fourier transform of the annihilation operators of a particular /-orbital. 



.7 = 1 



The operator in Eq. (38) satisfies the property 



'L 

IlMl, k=l. 



Thus a total momentum shift operator can be constructed as 

3=iin^^'^ (42) 

where 

2J(') =2J^'^^2J^')^. (43) 
The total momentum shift operator 3 can be used to construct a total position operator 



L-l 
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The operator Xpam includes the positions of electrons in impurity orbitals as well as those in the conduction band. To 
prove that it satisfies the three required conditions proceeds as before. Proving that the time-derivative of the position 
operator is equal to the current is simplified by the fact that the operators QT" commute with the periodic Anderson 
Hamiltonian. This is another consequence of the fact that there arc no hoppings between /-orbitals positioned on 
different sites. Hence all that needs to be proven is that the commutator corresponding to il gives the current operator 
corresponding to that of the Hubbard model [15]. This was already shown above. 

The operator X is well defined in the occupation number representation and it and its moments can thus be 
calculated in practical situations. Here we demonstrate the utility of the operator X by calculating the moments 
and performing finite size scaling for the Gutzwiller approximate solution of the Hubbard model at half-filling. The 
Gutzwiller wavefunction (GWF) has the form 

I*) = exp ^-T^^n^tnaj |*o). (45) 

where |\l/o) is a noninteracting wavefunction, and 7 is a variational parameter which projects out double occupations. 
Most often ]^'o) is the Fermi sea. In this case the exact solution in one [16,17] and infinite dimensions [18,19] are avail- 
able. At half- filling the former is metallic for finite U, in contradiction with the exact solution [20] . An approximate 
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solution to the GWF due to Gutzwillcr (GA) results in the Brinkman-Ricc metal insulator transition [11,21.22]. In fi- 
nite dimensions the GA is only approximate, however in infinite dimensions it correponds to the exact solution [18,19]. 
In a one- dimensional system the Brinkman-Rice transition is known to occur at Uc ~ 10. If \^o) is a non-interacting 
antifcrromagnctic wavcfunction the Gutzwillcr wavcfunction can be made insulating [23]. In the following, to assess 
the localization accompanying the metal-insulator transition wc calculate the quantity 



via quantum Monte Carlo methods [24,25]. 
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FIG. 1. X4 (defined in Eq. (46)) for the Hubbard model using the Gutzwiller wavefunction evaluated in the Gutzwiller 
approximation scheme. The Brinkman-Rice transition is known to occur at Uc ?s 10. 

In Fig. 1 X4 as a function of the Hubbard interaction strength for three different system sizes is presented. A 
transition at Uc ~ 10 is clearly visible from the simultaneous drop of all three curves. For large U {U > 11) the 
largest (smallest) system shows the smallest (largest) value of the fourth moment, which is the tendency one expects 
for the insulating state. (The same behaviour was found for the square-root of the second order deviation.) These 
results coincide with what is known about the Brinkman-Rice transition being a localization transition [22]. 




FIG. 2. Size dependence of X4 for a metallic Gutzwiller wavefunction. The inset shows the size dependence of the fourth 

order Binder cumulant. 

In Figs. 2 and 3 a metallic and an insulating wavefunction are compared. For the former the noninteracting 
wavefunctions (ground state of the U = Q system) was used in place of in Eq. (45). For the insulating 
wavcfunction an antiferromagnetic solution was used with a magnetization of m = 0.33333. The size dependence of 
the quantity Xi is clearly sensitive to whether the system is metallic or insulating: as the variational parameter 7 is 
increased Xi decreases in both cases, but the size dependence of X4 is opposite between the two cases. The metallic 
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state (Fig. 2) shows an increase in delocalization with system size, whereas in the insulating state (Fig. 3) the larger 
system is more localized. The insets in Figs. 2 and 3 show the value of the fourth order Binder cumulant [26-28] 
defined as 



a quantity used in the finite size scaling [29] of phase transitions. approaches a value of two-thirds in the case of 
perfect localization. Again, total order (localization) is approached by both the metallic and insulating wavefunctions, 
but the size dependence is the opposite between the two cases, with the larger system closer to the limiting value of 
two-thirds for the insulating wavefunction (hence more localized). 




FIG. 3. Size dependence of X4 for an insulating Gutzwiller wavefunction. The inset shows the size dependence of the fourth 

order Binder cumulant. 

In this paper a total position operator was derived for lattice models. The operator satisfies three crucial criteria: 
it is the generator of total momentum shifts, it is Hcrmitian, and its time derivative corresponds to the total current 
operator. The form of the operator is such that the average total position and its moments can bo readily calculated. 
Hence Binder ciimulants used in finite size scaling can also be evaluated. The sensitivity of such moments and 
cumulants was also demonstrated by investigating their size dependence in the Brinkman-Rice transition, and metallic 
and insulating Gutzwiller wavefunctions. 

Part of this work was performed at the Institut fiir Theoretische Physik at TU-Graz under FWF (Forderung der 
wissenschaftlichen Forschung) grant number P21240-N16. Part of this work was performed under the HPC-EUR0PA2 
project (project number 228398). Helpful discussions with H. G. Evertz are gratefully acknowledged. 



[1] W. Kohn, Phys. Rev., 133 A171 (1964). 

[2] E. I. Blount, Solid State Physics: Advances in Research and Applications, Eds. F. Seitz and D. Turnbull, 13 305 (1962). 

[3] N. D. M. Hine and W. M. C. Foulkcs, J. Phys.: Condens. matter, 19 506212 (2007). 

[4] R. Resta, Phys. Rev. Lett, 80 1800 (1998). 

[5] R. Rcsta and S. SorcUa, Phys. Rev. Lett, 82 370 (1999). 

[6] A. A. Aligia and G. Ortiz, Phys. Rev. Lett, 82 2560 (1999). 

[7] I. Souza, T. Wilkcns, and R. M. Martin, Phys. Rev. B, 62 1666 (2000). 

[8] J. Hubbard, Proc. Roy. Soc, A276 238 (1963). 

[9] J. Kanamori, Prog. Theoret Phys., 30 275 (1963). 
[10] M. C. Gutzwiller, Phys. Rev. Lett., 20 1445 (1963). 
[11] M. C. Gutzwiller, Phys. Rev., 137 A1726 (1965). 

[12] F. H. L. Essler, H. Frahm, F. Gohmann, A. Kliimper, and V. E. Korepin, The One-Dimensional Hubbard Model, Cambridge 

University Press, (2005). 
[13] G. D. Mahan, Many-Particle Physics, 3rd Ed., Kluwer Academic (2000). 



8 



[14] M. Imada, A Fujimori, and Y. Tokura, Rev. Mod. Phys., 70 1039 (1998). 

[15] D. Baeriswyl, C. Gros, and T. M. Rice, Phys. Rev. B, 35 8391 (1987). 

[16] W. Metzner and D. Vollhardt, Phys. Rev. Lett., 59 121 (1987). 

[17] W. Metzner and D. Vollhardt, Phys. Rev. B, 37 7382 (1988). 

[18] W. Metzner and D. Vollhardt, Phys. Rev. Lett, 62 324 (1989). 

[19] W. Metzner and D. Vollhardt, Helv. Phys. Acta, 63 364 (1990). 

[20] E. H. Lieb and F. Y. Wu, Phys. Rev. Lett, 20 1445 (1968). 

[21] W. F. Brinkman and T. M. Rice, Phys. Rev. B, 2 4302 (1970). 

[22] D. Vollhardt, Rev. Mod. Phys., 56 99 (1984). 

[23] W. Metzner, Z. Phys. B, 77 253 (1989). 

[24] H. Yokoyama and H. Shiba, J. Phys. Soc. Japan, 56 1490 (1986). 

[25] B. Hetenyi, H. G. Evertz, and W. von der Linden, Phys. Rev. B 045107 (2009). 

[26] K. Binder, Phys. Rev. Lett, 47 693 (1981). 

[27] K. Binder, Ferroelectrics, 73 43 (1987). 

[28] K. Binder, Annu. Rev. Phys. Chem., 43 33 (1992). 

[29] M. E. Fisher and M. N. Barber, Phys. Rev. Lett, 28 1516 (1972). 



9 



